#! /bin/sh
# Velocity analyses for the cmp gathers
# Authors: Dave Hale, Jack K. Cohen, with modifications by John Stockwell
# NOTE: Comment lines preceeding user input start with  ##
#set -x

## Set parameters
datadir=.

velpanel=gain.spike.mute.pbal.cdp.su   # gained data sorted in cdps
vpicks=radonnmovel.par       # output file of vnmo= and tnmo= values
normpow=0               # see selfdoc for suvelan
slowness=0              # see selfdoc for suvelan
cdpfirst=1              # minimum cdp value in data
cdplast=2142            # maximum cdp value in data
cdpmin=128              # minimum cdp value used in velocity analysis
cdpmax=2110             # maximum cdp value used in velocity analysis
dcdp=512                # change in cdp for velocity scans
fold=120                # maximum number of traces per cdp gather
dxcdp=12.5              # distance between successive midpoints
                        # in full data set
mix=0			# number of adjacent cdp panels on either side
			# of a given panel to mix

## Set velocity sampling and band pass filters
nv=200			# number of velocities in scan
dv=15.0			# velocity sampling interval in scan	
fv=1450.0		# first velocity in scan
nout=1500		# ns in data		

## Set interpolation type for velocity function plots
interpolation=mono	# choices are linear, spline, akima, mono

## set filter values for wiggle trace plots
f=1,5,70,80		# bandwidth of data to pass
amps=0,1,1,0		# don't change

## suximage information
wclip=0		# This number should be between 0 to .15 for real data 
bclip=.3	# this number should be around .5
cmap=hsv2	# colormap 
perc=97		# clip above perc percential in amplitude
xcur=2		# allow xcur trace xcursion in wiggle trace plots
curvecolor=black	# color of stacking velocity picks curve


#average velocity
vaverage=2100        # this may be adjusted

# radon transform parameters
dp=8				# increment in p in radon transform
pmin=-2000			# minimum value of p in radon transform
pmax=2000			# maximum value of p in radon transform
pmula=90			# t=tmax intercept of radon filter
pmulb=90				# t=0 intercept of radon filter
offref=-3237			# offset at maximum moveout
interoff=-262			# offset at minimum moveout

######## You shouldn't have to change anything below this line ###########

# binary files output
vrmst=vrmst.bin		# VRMS(t) interpolated rms velocities
vintt=vintt.bin		# VINT(t,x) as picked
vinttav=vinttav.bin	# average  VINT(t) of VINT(t,x)
vrmstav=vrmstav.bin	# average  VRMS(t) of VRMS(t,x)
vinttuni=vinttuni.bin	# interploated Vint(t,x)
vintzx=vintzx.bin	# VINT(z,x)interpolated interval velocities
vintzav=vintzav.bin	# average  VINT(z) of VINT(z,x)
vintxz=vintxz.bin	# VINT(x,z)interpolated interval velocities

### Get header info
cdpcount=0		 #  counting variable
dxout=0.004		# don't change this

nout=`sugethw ns <$velpanel | sed 1q | sed 's/.*ns=//'`
dt=`sugethw dt <$velpanel | sed 1q | sed 's/.*dt=//'`
dxout=`bc -l <<END
	$dt / 1000000
END`

cdptotal=`bc -l <<END
	$cdplast - $cdpfirst
END`

dtsec=`bc -l <<END
        $dt / 1000000
END`

echo  "Skip Introduction? (y/n) " | tr -d "\012" >/dev/tty
read response
case $response in
n*) # continue velocity analysis


### give instructions
echo
echo
echo "            Instructions for Velocity Analysis."
echo
echo " This version of the Velan script is expressly for picking velocities"
echo " to use in radon transform domain multiple suppression. Picks are"
echo " to be made below and to the right of the maxima that are normally"
echo " picked for NMO-Stack."
echo
echo " An  image semblance map will appear on the right of your screen."
echo " A wiggle trace plot of the cdp panel being analysed will appear"
echo " on the left as an aid in picking. Click on the semblance contour"
echo " map to make that window active."
echo
echo " Pick velocities by placing the cursor on each peak (do not click) "
echo " in the semblance plot and typing 's'. Type 'q' when last peak is "
echo " picked. Picked values will be saved into a file."
echo
echo " Note, this is 'blind' picking. You will not see any indication"
echo " on the contour plot that a point has been picked."
echo
echo " Note also, that it is best if a value of the velocity is picked "
echo " at the beginning of the data (t=0.0 usually). The picks must "
echo " be increasing in time.  If you feel you have made an incorrect pick"
echo " you will be given an opportunity to pick the velocities again. "

pause

echo
echo
echo " You will be shown a wiggle trace plot and a semblance map of the"
echo " multiple-suppressed data. Do not attempt to pick these. You will"
echo " be asked if the multiple-suppression is ok. If you will answer y"
echo " you will be asked if you want to continue doing velocity picking."
echo " If you answer y, the script will take you to the next cdp."
echo " If you answer n, you will be asked if you want to perform picking"
echo " on the multiple suppressed data. If you answer n you the picking"
echo " will be purformed on the non-multiple suppressed panel."
echo
echo " The radon processing is always applied on the original panel, so"
echo " you can back off on your velocities, if necessary, and the next"
echo " set of plots will reflect this."

echo
echo " A graph of the velocity function will appear, and a prompt to" 
echo " hit the return key will be seen in this terminal window.  You"
echo " will then see an nmo corrected version of the cdp gather you that"
echo " you are performing velocity analysis on." 
echo
echo " Caveat, the interval velocities plotted should not be taken too"
echo " seriously as these are interpolated from just a handful of values."
echo " This curve is supplied for reference only."
echo
echo " You will be prompted in the terminal window to hit return. Then "
echo " you will be  will be asked if your picks are ok. This gives you "
echo " a chance to re-pick the velocities if you do not like the velocity"
echo " function you have obtained."
echo
echo " When the shell script reruns, it will ask you about reusing velocity"
echo " panels. Once these are made, you do not need to remake them, so answer"
echo " yes to reuse, in this case."
echo
echo
echo " Advanced features"
echo
echo " You may set the mix=  to a value that will indicate the number of"
echo " cdps to the right and left of the current cdp to \"mix\" or average"
echo " with the current cdp on which velocity analysis is being performed".
echo

;;
*y) #continue

echo
echo
echo "Beginning the velocity analysis"
echo
echo
echo

;;
esac

########################### start velocity analysis #####################


ismix=false
cdp=$cdpmin
while [ $cdp -le $cdpmax ]
do
	cdpcount=` expr $cdpcount + 1 `
	ok=false
	reusepanel=false

	
	# see if panel.$cdp exists
	if [ -f $datadir/panel.$cdp ]
	then
		echo  "panel.$cdp exists. Reuse? (y/n) " | tr -d "\012" >/dev/tty
		cp $datadir/panel.$cdp $datadir/origpanel.$cdp
		read response
		case $response in
		n*) # continue velocity analysis
			reusepanel=false
		;;
		y*) # no need to get velocity panel
			reusepanel=true
		;;
		esac
	fi

	# see if par.$cdp and $vrmst.$cdp exist
	if [ -f par.$cdp ]
	then

		if [ -f $vrmst.$cdp ]
		then
			echo
			echo " file $vrmst.$cdp already exists"
			echo " indicating that cdp $cdp has been picked"
		fi
		echo
		echo " file par.$cdp already exists"
		echo " indicating that cdp $cdp has been picked"
		echo

		echo  "Redo velocity analysis on cdp $cdp? (y/n) " | tr -d "\012" >/dev/tty
		read response
		case $response in
		n*) # continue velocity analysis with next cdp
			ok=true
		;;
		y*) # continue with same value of cdp
			ok=false
		;;
		esac
	fi

	# begin velocity analysis
	while [ $ok = false ]
	do
		echo "Starting velocity analysis for cdp $cdp"

		touch npairs.$cdp curvefile.$cdp

		# capture the current cmp gather
		if [ $reusepanel = false ]
		then
			rm temp temp1
			if [ $mix -ne 0 ]
			then
				cdplower=` expr $cdp - $mix `
				cdpupper=` expr $cdp + $mix `

				echo $mix $cdp $cdplower $cdpupper

			        suwind < $velpanel key=cdp \
					min=$cdplower \
					max=$cdpupper > $datadir/mixpanel.$cdp

				susort offset cdp < $datadir/mixpanel.$cdp > $datadir/temp1
				sustack key=offset < $datadir/temp1 > $datadir/temp 	
				sushw key=cdp a=$cdp < $datadir/temp > $datadir/panel.$cdp
				cp $datadir/panel.$cdp $datadir/origpanel.$cdp
				reusepanel=true
				ismix=true
			elif [ $mix -eq 0 ]
			then
				suwind < $velpanel key=cdp min=$cdp max=$cdp \
					count=$fold > $datadir/panel.$cdp 
				cp $datadir/panel.$cdp $datadir/origpanel.$cdp
				reusepanel=true
			fi
		fi

		# plot the wiggle traces of the cmp gather
		suxwigb < $datadir/panel.$cdp  title="CDP gather for cdp=$cdp" \
			ybox=10	xbox=50  key=offset \
				perc=$perc xcur=$xcur wbox=300 &

		# make a semblance panel
		sufilter < $datadir/panel.$cdp f=$f amps=$amps |
		suvelan nv=$nv dv=$dv fv=$fv |
		suximage wclip=$wclip bclip=$bclip \
		par=npairs.$cdp curve=curvefile.$cdp curvecolor=$curvecolor \
		f2=$fv d2=$dv ybox=10 xbox=450 wbox=600 legend=1 \
		units="semblance" cmap=$cmap \
		label1="Time (sec)" label2="Velocity (m/sec)" \
		title="Velocity Scan (semblance plot) for CMP $cdp" \
		mpicks=mpicks.$cdp

		# capture curve file
		cp mpicks.$cdp curvefile.$cdp
		a2b n1=2 outpar=junk < mpicks.$cdp > mpicks_bin.$cdp
		sed 's/n2=/npair=/g' < junk > npairs.$cdp


		# make a par file for this cdp based on the velocity and time picks
		sort <mpicks.$cdp  -n |
		mkparfile string1=tnmo string2=vnmo >par.$cdp

		# view the picked velocity function  as both inter
		echo "Putting up velocity function for cdp $cdp"
		sed <par.$cdp '
			s/tnmo/xin/
			s/vnmo/yin/
		' >unisam.p
		unisam nout=$nout fxout=0.0 dxout=$dxout \
			par=unisam.p method=$interpolation > uni.temp.bin
		velconv intype=vrmst outtype=vintt nt=$nout dt=$dt < uni.temp.bin > vel.temp.bin
		cat uni.temp.bin vel.temp.bin > graph.temp.bin

		zap xwigb 
		zap ximage 

		pause 
		echo
		echo
		echo
		echo
		echo
		echo "Wait a moment while radon is applied"
		echo "do not hit any keys until you are prompted."
		echo
		echo
		echo
		echo
		echo

		# view a semblance panel of the data filtered in the
		# radon domain
		## nmo->radon-> inverse NMO: for multiple suppression

		# mixing adjacent cdps?
		if [ $ismix = true ]
		then
			cp $datadir/mixpanel.$cdp  $datadir/temp
		elif [ $ismix = false ]
		then
			cp $datadir/origpanel.$cdp $datadir/temp
		fi

			
		# apply radon transform 
		sufilter < $datadir/temp f=$f amps=$amps |
		sunmo par=par.$cdp  |
		suradon offref=$offref interoff=$interoff pmin=$pmin \
		pmax=$pmax dp=$dp choose=1 igopt=2 \
		pmula=$pmula pmulb=$pmulb \
		depthref=1000 | sunmo par=par.$cdp invert=1 > $datadir/radonpanel.$cdp

		# perform mix and stack
		if [ $ismix = true ]
		then
			cp $datadir/radonpanel.$cdp $datadir/temp
			susort offset cdp < $datadir/temp > $datadir/temp1
			sustack key=offset < $datadir/temp1 |
			sushw key=cdp a=$cdp > $datadir/radonpanel.$cdp
		fi
		
		# show velan panel result
		suvelan < $datadir/radonpanel.$cdp nv=$nv dv=$dv fv=$fv |
		suximage wclip=$wclip bclip=$bclip \
		par=npairs.$cdp curve=curvefile.$cdp curvecolor=$curvecolor \
		f2=$fv d2=$dv xbox=620 wbox=600 legend=1 \
		units="semblance" cmap=$cmap \
		label1="Time (sec)" label2="Velocity (m/sec)" \
		title="DO NOT PICK! Test velocity scan after multiple suppression $cdp" &
		
		# show original panel
		suxwigb < $datadir/origpanel.$cdp \
			 title="Original data cdp=$cdp" \
				xbox=10 ybox=10  key=offset \
				perc=$perc xcur=$xcur wbox=300 &

		# show radon multiple suppressed panel
		suxwigb < $datadir/radonpanel.$cdp \
			 title="Multiple-suppressed cdp=$cdp" \
				xbox=320 ybox=10 key=offset \
				perc=$perc xcur=$xcur wbox=300 &

		echo
		echo
		echo
		echo  "View, do not pick!"
		echo
		echo
		echo

		echo
		echo
		echo
		pause
		echo
		zap xwigb
		echo
		echo
		echo
		echo "           Velocities and NMO corrected panel "

		xgraph < graph.temp.bin n=$nout nplot=2 d1=$dxout \
			linewidth=2,3  f1=0.0 width=250 height=700 \
			label1="Time (sec)" label2="Velocity (m/sec)" \
			title="Red:RMS Velocity function  Blue: approx Interval Velocity: CMP $cdp" \
			linecolor=2,4  &

		# view an  the multiple-suppressed panel no NMO
		echo "Hit return after nmo panel comes up"
		suxwigb < $datadir/radonpanel.$cdp \
			 title="Multiple-suppressed no nmo cdp=$cdp" \
			ybox=10	 key=offset \
				perc=$perc xcur=$xcur xbox=280 wbox=250 &

		# show NMO of panel with no additional multiple suppression
                sunmo par=par.$cdp < $datadir/panel.$cdp |
                suxwigb \
		title="NMO of cdp=$cdp (no multiple suppression)" \
		 xbox=500 wbox=250 xcur=$xcur \
		   ybox=10 perc=$perc xcur=$xcur  key=offset &

		# shown NMO with multiple suppression
                sunmo par=par.$cdp < $datadir/radonpanel.$cdp |
                suxwigb \
		title="NMO of cdp=$cdp (with multiple suppression)" \
		 xbox=750  wbox=250 xcur=$xcur \
		   ybox=10 perc=$perc xcur=$xcur  key=offset &

		pause  

		# kill windows
		zap xgraph 
		zap xwigb
		zap ximage

		echo
		echo
		echo
		# check to see if the picks are ok
		echo "If this is the first iteration of velocity picking,"
		echo " answer \"n\" to the next question. "
		echo  "Suppression and velocity picks  OK? (y/n) " | tr -d "\012" >/dev/tty
		read response
		case $response in
		n*) ok=false 
			echo
			echo
			echo
			echo "Picking again"
			echo
			echo
			echo
			# move the radon filtered panel and repeat
			# velocity analysis on that panel
			echo "Repeat analysis, this time using the multiple-suppressed panel? (y/n)"|
				tr -d "\012" >/dev/tty
			read response
			case $response in
			n*) # get a fresh copy of panel.$cdp
				# do nothing
				zap ximage
				zap xwigb

				cp $datadir/panel.$cdp origpanel.$cdp
                        	suwind < $velpanel key=cdp min=$cdp max=$cdp \
                                        count=$fold > $datadir/panel.$cdp

			;;
			*) # move radon panel into cdp panel
				mv radonpanel.$cdp panel.$cdp
				zap ximage
				zap xwigb
			;;
			esac
			
		;;
		*) ok=true 
			# capture resampled velocity
			unisam nout=$nout fxout=0.0 dxout=$dxout \
			par=unisam.p method=$interpolation > $vrmst.$cdp

			# clean up the screen
			zap ximage
			zap xgraph
			zap xwigb
			zap xcontour
		;;
		esac

	done </dev/tty

	echo
	echo
	echo  "Continue with velocity analysis? (y/n) "  | tr -d "\012" >/dev/tty
	read response
	case $response in
	n*)	# if quitting set cdp to a value large enough to
		# break out of loop 
		cdp=`expr $cdpmax + 1`
	;;
	y*)
		# else get the next cdp
	cp curvefile.$cdp tempcurve
	cp npairs.$cdp tempnpairs
	cdp=`bc -l <<END
		$cdp + $dcdp
END`
	mv tempcurve curvefile.$cdp
	mv tempnpairs npairs.$cdp

	;;
	esac

done

set +x


### Combine the individual picks into a composite sunmo par file
echo "Editing pick files ..."
>$vpicks
echo  "cdp=" | tr -d "\012" >>$vpicks
cdp=$cdpmin
echo  "$cdp"  | tr -d "\012" >>$vpicks
cdp=`bc -l <<END
	$cdp + $dcdp
END`
while [ $cdp -le $cdpmax ]
do
	echo  ",$cdp"  | tr -d "\012" >>$vpicks
	cdp=`bc -l <<END
		$cdp + $dcdp
END`
done
echo >>$vpicks

cdpcount=0
rm $vrmst
cdp=$cdpmin
while [ $cdp -le $cdpmax ]
do
	cat $vrmst.$cdp >> $vrmst
	cat par.$cdp >>$vpicks
	cdp=`bc -l <<END
		$cdp + $dcdp
END`
	cdpcount=` expr $cdpcount + 1 `
done

# build velocity files to be used for comparison (not recommended for migration)
vrmstpar=vrmst.par
vinttpar=vintt.par
vinttplotpar=vinttplot.par
vintzplotpar=vintzplot.par
unipar=unisam.par

# build par files
echo "n1=$nout d1=$dt n2=$cdpcount f2=$cdpmin d2=$dcdp " > $vrmstpar
echo "nt=$nout d1=$dt ns=$nout nx=$cdpcount fx=$cdpmin dx=$dcdp  " > $vinttpar
echo "n=$nout nplot=1 d1=$dxout x2beg=$fv style=seismic width=400 height=700  " > $vinttplotpar
echo "nx1=$nout nx2=$cdpcount n1=$nout n2=$cdptotal" > $unipar

# convert rms velocities to interval velocities 
velconv intype=vrmst outtype=vintt  par=$vinttpar < $vrmst > $vintt

# make an average velocity profile
suaddhead < $vintt ns=$nout | sustack | sustrip  > $vinttav 

# make an average velocity profile
suaddhead < $vrmst ns=$nout | sustack | sustrip  > $vrmstav 

# build a uniformly sampled v(t,x) velocity profile
unisam2 < $vintt par=$unipar  | smooth2 r1=5 r2=5 par=$unipar >  $vinttuni

                                                                                
# get depth sampling interval
dzout=`bc -l <<END
        ( $vaverage * $dtsec ) / 2.0
END`
                                                                                
echo $dzout

# par file for average interval velocity
echo "n=$nout nplot=1 x2beg=$fv d1=$dzout style=seismic xbox=500 width=400 height=700  " > $vintzplotpar
                                                                                
# build v(z,x) 
velconv intype=vintt outtype=vintz  dt=$dtsec \
nx=$cdplast nz=$nout dz=$dzout < $vinttuni |
smooth2 r1=10 r2=20 n1=$nout n2=$cdplast > $vintzx

suaddhead < $vintzx ns=1500 |
sushw key=cdp a=1 |
sustack | sustrip > $vintzav
                                                                                
# build v(x,z)
transp < $vintzx n1=$nout > $vintxz



# final echos
echo "V(t) RMS (stacking) velocity file: $vrmst is ready"
echo "V(t,x) Interval velocity file: $vintt is ready"
echo "V(z) INT (interval) velocity file: $vintz is ready"
echo "V(z,x) Interval velocity file: $vintzx is ready"
echo "V(x,z) Interval velocity file: $vintxz is ready"
echo "sunmo par file: $vpicks is ready"


exit 0
